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The equilibrium ensemble approach to disordered systems is used to investigate the critical be- 
haviour of the two dimensional Ising model in presence of quenched random site dilution. The 
numerical transfer matrix technique in semi-infinite strips of finite width, together with phenomeno- 
logical renormalization and conformal invariance, is particularly suited to put the equilibrium en- 
semble approach to work. A new method to extract with great precision the critical temperature 
of the model is proposed and applied. A more systematic finite-size scaling analysis than in pre- 
vious numerical studies has been performed. A parallel investigation, along the lines of the two 
main scenarios currently under discussion, namely the logarithmic corrections scenario (with crit- 
ical exponents fixed in the Ising universality class) versus the weak universality scenario (critical 
exponents varying with the degree of disorder), is carried out. In interpreting our data, maximum 
care is constantly taken to be open in both directions. A critical discussion shows that, still, an 
unambiguous discrimination between the two scenarios is not possible on the basis of the available 
finite size data. 
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I. INTRODUCTION 

Recent years have witnessed renewed efforts in the sta- 
tistical physics community to understand phase transi- 
tions in simple disordered classical spin systems. As a 
matter of fact, these efforts have produced conflicting 
statements concerning the effects of disorder on critical 
phenomena, covering almost the complete spectrum of 
conceivable alternatives. 

This holds in particular, but not exclusively, for two- 
dimensional {2d) disordered ferromagnetic (i.e., unfrus- 
trated) Ising models. These models have been widely 
studied, both because the corresponding pure system is 
well understood and because they constitute a marginal 
case in the Harris criterion^ which assesses whether disor- 
der constitutes a relevant or irrelevant perturbation for 
the critical behaviour of the pure system. For models 
of this type, the discussion currently appears to nar- 
row down on two conflicting scenarios, namely the loga- 
rithmic corrections^^^'^ versus the weak universality^^~'^^ 
scenario, though a broader spectrum of alternatives had 
been discussed earlier^^^^^ (the interested reader will find 
a comprehensive report on the literature up to approxi- 
mately 1982 in an early review by Stinchcombe^^ ) . We 
will describe these scenarios and discuss these (and re- 
lated) results in greater detail later on. 

The object of our study is the randomly spin di- 
luted 2d Ising model. Our investigation is based on two 
main ingredients. First, we use the equilibrium ensem- 
ble approach to disordered systems^^'^^'^*'^^'^* to map 
the quenched system onto an equivalent thermodynamic 
equilibrium system in an enlarged phase space. Perform- 
ing this mapping exactly would be tantamount to pro- 
viding an exact solution to the original problem, which 
is clearly infeasible. So a scheme of approximations based 
on a moment matching idea is invoked. Second, we re- 
sort to conventional transfer matrix (TM) techniques to 
implement the method in practice on finite-width strips, 
analyzing finite-size results in the line of Nightingale's 



phenomenological renormalization group scheme. ^^'^^ 

The purpose of the present paper is to provide details 
of our TM study,^*'^^ as well as to include new material 
and to describe significant advances in the understanding 
of finite-size scaling (FSS) signatures in the presence of 
logarithmic corrections, which together have allowed us 
to boost the accuracy of our results considerably and to 
obtain a sounder appreciation of the subtleties that may 
emerge in the interpretation of the data. 

The main and unexpected finding in Ref. 21 has been 
a continuous variation of the critical exponents a, /?, 7 
and v with the spin density p in a manner which was 
observed to comply with the idea of weak universality.^^ 
That is, the exponent rj describing the decay of critical 
correlations, and the magnetic exponent 6, as well as the 
ratios /3/i^ and 'y/iy were found to be independent of p. 
The results were obtained by extrapolations of FSS data 
based on rather moderate strip-widths, and they were 
in complete quantitative and qualitative agreement with 
those of a Monte Carlo study by Kim and Patrascioiu.^" 
These results are in conflict with those supporting the 
logarithmic corrections scenario, where modifications of 
the relevant thermodynamical quantities at criticality ap- 
pear through logarithmic terms in the reduced temper- 
ature, while the system is left in the same universality 
class (i.e. with the same critical exponents) as that of 
the pure 2d Ising model. These findings are almost all 
concerned with the ferromagnetic bond disordered situa- 
tion, either theoretically, '^"^'^^ or via numerical^'^'^"-'^^'-'^^ 
and series expansion^^ approaches, but the comparison 
with the spin diluted model is possible supposing ~ as it 
is generally believed to be the case - that the two systems 
are in the same universality class. ^^'^^ 

In discussing our results, we have to address two main 
issues. The first is concerned with the reliability of our 
method, which is based on an approximate description 
of quenched disorder. The evidence we have been able 
to compile does give us strong confidence in the validity 
of our approach. This granted, we turn to the second 
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issue: can - or more precisely to what extent can - our 
results provide evidence in favour of or against any of 
the conflicting scenarios so far advanced to describe the 
critical behaviour of 2d disordered ferromagnetic Ising 
models? This is indeed a subtle question, and we devote 
almost two Sections (Sees. IV and V) to discussing it. 
It turns out that a considerable portion of the available 
finite-size data from simulations or TM studies - includ- 
ing those presented in our earlier study^^'^^ as well as 
some new ones - may perhaps not allow to decide with 
sufficient confidence between the two most serious can- 
didates. And we shall explain why. Larger system sizes 
would in any case be necessary in studies along those 
lines to permit taking sides. 

Finally, since one of the aims of our study is to dis- 
criminate between these two contradicting scenarios, one 
of our main concerns is an open-minded attitude, per- 
haps not easily found in the literature, in the analysis of 
the numerical results: we take the greatest care not to 
select one of them a-priori and to fit our data to it, but 
we rather attempt to be open in both directions. 

The remainder of our paper is organized as follows. In 
Sees. II A and II B, in order to set the scene and to fix our 
notation, we describe the model that is studied in the se- 
quel and we briefly recall our methodical background, viz. 
the equilibrium ensemble approach to disordered systems 
(EEA). While Sec. II C presents some details related to 
the computational side of the problem, we reserve Sec. 
Ill for a brief review of the two above mentioned sce- 
narios, and Sec. IV for a presentation and discussion of 
our main results. We describe a new finite-size proce- 
dure to extract the critical temperature of the model as 
a function of the degree of dilution, which enables us to 
draw an accurate approximation of the phase diagram. 
We then determine values for the critical exponents, for 
the central charge and for the specific heat of our ap- 
proximating systems. We end our paper with a critical 
report of the available finite-size data obtained also by 
other authors, and with a comprehensive discussion both 
on the method and on the outcome of our investigation 
(Sec. V). An Appendix, mainly based on renormaliza- 
tion group results by Cardy and Ludwig,'^^^'''' presents a 
systematic FSS analysis of thermodynamic quantities in 
the presence of logarithmic corrections. 

II. SETTING THE SCENE 

A. The system under investigation 

The object of our investigation is the 2d randomly spin 
diluted Ising model, described by the Hamiltonian 

H{a\K) = —J^^kiaikjaj — H^^^i^^i ■ (1) 

{ij) » 

Here, the first sum is over all nearest neighbour pairs 
(ij) of, say, a square lattice A, and the second over all 
sites i. The ai denote Ising spins and the ki G {0, 1} 
are occupation numbers signifying whether in a disorder 



configuration k a site i is occupied by a spin {ki = 1) 
or not {ki = 0). The ki are taken to be quenched ran- 
dom variables, i.e., they are fixed, and randomly chosen 
according to the probability 

g(/c)=n/'(i-p)'-'' , (2) 

i 

which simply requires each site to be occupied with prob- 
ability yO, and to be empty with probability 1 — p. Thus 
p is the average density of spins in the system. 

Our aim is to study the thermodynamics of the system, 
described by the (dimensionless) quenched free energy 

f, = -N-^{\nZN{K))^ , (3) 

i.e., the average of the system's free energy over the distri- 
bution q{K) describing the statistics of the disorder con- 
figurations K. Here TV denotes the system size, and Zn{k) 
is the partition function of a system of size |A| = A'' at 
fixed disorder configuration k. In particular, we would 
like to locate the phase-transition line Tc{p) below which 
the system is known to develop spontaneous ferromag- 
netic order in the thermodynamic limit N ^ oo, and to 
study the critical behaviour of the system at and around 
this line. Whereas the thermodynamics of the Id system 
is known exactly, an exact solution of the 2d model is cur- 
rently way beyond our abilities. Approximation methods, 
which may, however, well produce exact results in certain 
limiting cases are thus called for. 

B. Methodical background: the equihbrium 
ensemble approach 

The EEA to disordered systems, as well as its appli- 
cation to the 2d spin diluted Ising model, have been ac- 
curately described elsewhere. Therefore we will not 
enter again into the details of the method, rather just 
concisely recall the main ideas. 

The central idea is to treat the configurational degrees 
of freedom ki in a system with quenched disorder on the 
same footing as the dynamical variables proper - the 
spins, in the case at hand. That is, the phase space 
is enlarged to include the occupation-number configura- 
tions K, and an additional term depending only on k, the 
so-called disorder potential (f){K) is added to the Hamil- 
tonian 

H{a\K) n'^{(j, k) ^ H{a\K) + 0(k) , (4) 

which is then determined in such a way that configuration 
averaging, as implied by Eq. (3), becomes part of the 
Gibbs average in the enlarged phase space. 

To utilize the method in practice, one still needs an 
explicit representation for 4>{k). For spin or site diluted 
systems the following representation in terms of products 
of occupation numbers 

-l3(j){K) = Ai Y^{h - p) + Aa Y^{hkj -p^)+ ... 

i {ij) 
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is always suitable. The first sum explicitly displayed in 
(5) is over all lattice sites, the second over all nearest 
neighbour pairs, the third over all elementary plaquettes 
(of size |P|) of the system, and so on. 

It has been realized^^ that the self-consistent con- 
straint equations which determine the couplings \v as a 
function of temperature T, magnetic field H and spin 
density p (Eqs. 4 of Ref. 21), can be written as nec- 
essary conditions describing the maxim,um of the equi- 
librium ensemble's dimensionless free energy (per site) 
= -N-HnZ'l' = -N-HnJ2a^^exp[-j3n'l'{a,K)] 
with respect to the X,/s. They can therefore be written, 
and are in practice evaluated too, as 



df-P_ 
' 3X2 



ih)^ - P = , 



dX, 



pl^l=0, ... . 



(6) 



Solving Eqs. (6) for every cluster of connected sites 
- i.e. exactly obtaining the quenched free energy via 
the EEA - is clearly infeasible. However, a systematic 
scheme of approximations may be put up by matching 
only a subset of the full set of moments of the equilib- 
rium c;nsemble's Gibbs distribution to the corresponding 
subset of moments of the problem with quenched disor- 
der. This amounts to forcing ki correlations on larger and 
larger groups of lattice sites to coincide with those of the 
quenched system. The larger the number of constraints 
taken into account, the better a description of the fully 
quenched system is obtained. 

In Ref. 21 four different approximating systems named 
(a) - (d) were studied, differing by the number of terms 
kept in (5), and thus by the sot of constraints which is 
actually retained in (6). The present study is focused 
only on approximation (d): the first three contributions 
to the disorder potential are kept, i.e. single site, pair 
and square plaquctte terms. In addition, due to the 
anisotropy of the strip geometry employed in the calcu- 
lation (see Sec. II C below) , pair occupancy parallel and 
perpendicular to the strip are treated as separate con- 
straints, i.e. they arc controlled by separate couplings 
A211 and A2_L in the disorder potential which, in turn, 
have to be determined by two distinct equations of the 
type (6). 



C. Transfer Matrix Approach 

The approach outlined above might at least in princi- 
ple - solve our problem to any desired degree of accuracy. 
However, in > 2 not even the simplest approximating 
system is exactly solvable. To analyse the thermody- 
namics or the critical behaviour of our systems, we use 
the strip FSS techniques pioneered by Nightingale.^^ A 
distinct advantage of the EEA in this context is that 
we are dealing with translationally invariant equilibrium 



systems with short range couplings at all levels of approx- 
imation within the moment matching scheme described 
above. FSS can therefore be implemented using conven- 
tional TM techniques for nonrandom systems. 

Our strip is a. L x L' square lattices, with the ther- 
modynamic limit taken in the L' direction, and periodic 
boundary conditions (unless otherwise stated) imposed 
in the perpendicular direction. The row to row transfer 
matrix F, constructed from the local Boltzmann factors, 
can be chosen to be symmetric for all systems (a) - (d) 
of Ref. 21, and is diagonalized for strips of fixed width L, 
with the appropriate set of constraints (Eqs. (6)). 

The dimensionless free energy (per site) is related 
to the largest eigenvalue 71 of F via 



-/L = L-^ln7i 



(7) 



(to simplify notation, we omit stating the L dependence 
of F and its eigenvalues in what follows). 

The constraint equations (6) have been formulated in 
terms of first derivatives of the free energy, and thereby 
in terms of first derivatives of f^. As mentioned above, 
those derivatives are nothing but gradient information in 
the problem of maximizing (or minimizing — /t) over 
the appropriate set of couplings. Notice that the tran- 
scendental nature of these equations necessitates an iter- 
ative numerical solution. Having in this way determined 
the (approximate) disorder potential, one may compute 
thermodynamic functions and, in particular, the corre- 
lation length of the order parameter fluctuations, which 
is needed in the phenomenological renormalization group 
scheme. The correlation length is given in terms of the 
largest and the second largest eigenvalue of the TM as 



In 



72 
71 



(8) 



Thermodynamic functions of interest are obtained by dif- 
ferentiating the free energy with respect to temperature 
or magnetic field. Here, care must be taken because of 
implicit field and temperature dependences in the cou- 
plings of the disorder potential. In other words, deriva- 
tives have to be performed along the solution-manifold of 
the appropriate set of constraint equations. 

Introducing modified spins Si = fcitXi with values in 
{0, ±1}, in a single variable we encode, both, presence or 
absence of a spin (the occupation numbers being given 
by ki = \Si\), and the spin states. Since the dimension 
of F grows as 3^ x 3^, we have used two different strate- 
gies - with complementary strengths and weaknesses - 
to simplify the computational task of dealing with huge 
transfer matrices. The first is a group theoretical analy- 
sis which exploits the various symmetries of the matrix 
and effectively reduces its dimensionality. The second is 
the use of sparse matrix techniques. 

We will not bother the reader by giving a detailed ac- 
count of these techniques since, apart from the fact that 
they have been non-trivially specialized to deal with our 
particular system, they are not new (for the sparse ma- 
trix technique see, e.g., Ref. 30). We just notice that 
within either of them it is possible to identify the blocks 
F^^^ and F^^', in the block diagonal representation of F, 
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that contain its two largest eigenvalues 71 and 72 re- 
spectively, r^^^ is the block spanned by the space of 
eigenvectors transforming symmetrically under spin re- 
versal, while 72 is found in the block spanned by the 
eigenvectors of T which are anti-symmetric under spin 
reversal. In a similar spirit, one may enquire into the 
large distance behaviour of the ki correlation functions 
Gk{rij) = {hkj) — {ki){kj). This requires to locate a cor- 
responding eigenvalue 72, which controls the asymptotic 
decay of Gk{r). It is found as the second largest eigen- 
value in the block corresponding to a representation of 
the symmetry group which is symmetric with respect to 
both external (spatial) and internal (spin) symmetries of 
the system. Indeed, whereas the spins Ui change sign 
under global spin reversal, entailing that the eigenvalue 
to be inserted in (8) for the computation of the spin- 
spin correlation length must belong to the antisymmetric 
block of r, the disorder variables fc, = l^jl do not. 

The computation of second order derivatives of 71, 
needed to evaluate some of the relevant thermodynam- 
ical quantities, is feasible only in the group theoretical 
reduction approach. The second derivatives in fact re- 
quire knowledge of the complete eigensystem, which is 
still far too large within the sparse matrix approach. In 
this case one has to resort to finite differences of first or- 
der derivatives to compute, e.g., the specific heat or the 
spin susceptibility. The advantage of the sparse-matrix 
approach, on the other hand, is that it can be pushed 
to considerably larger strip-widths than the group the- 
ory approach. We reach the value L = 9 with the latter 
and L = 13 with the former. Moreover, this approach 
does not exploit any symmetries beyond the spin rever- 
sal symmetry in zero magnetic field. It is therefore also 
applicable if one introduces an anti-ferromagnetic seam 
along the strip, a device which can be used to compute 
the domain wall free energy We shall exploit this quan- 
tity below to considerably boost the precision of our FSS 
analysis. 



III. THE TWO SCENARIOS 

Before showing our results, it is time to briefiy present 
the two main scenarios which have survived over a wider 
range of possibilities for the description of the critical 
behaviour of the 2d site diluted Ising model, but which 
are obviously mutually excluding. 

The logarithmic corrections scenario is based on the 
quantum field theory results by Dotsenko and Dotsenko^ 
and Shalaev,^ Shankar^ and Ludwig^ (the latter contri- 
butions correcting certain errors in the former). Though 
strictly valid only in the limit of weak disorder, they in- 
dicate that the presence of impurities affects the critical 
properties of the model only through a set of logarithmic 
corrections to the pure system behaviour. In particular, 
according to this picture, the correlation length of the 
infinite system close to the phase transition is expected 
to show the form 



where t is the reduced temperature t = {T — Tc)/Tc ^ 1 
(Tc being the critical temperature, which will depend on 
p), the exponents v and i> are respectively 1 and 1/2, and 
5 is a non- negative constant, such that ^ = in the pure 
case (usual power law behaviour recovered) and increases 
with increasing disorder. A similar behaviour holds for 
the magnetic susceptibility: 



Xo 



l + gln( i 



(10) 



with 7 = 7/4 and 7 = 7/8. For the specific heat, too, 
the exponent of the pure case a = is not modified by 
the introduction of disorder, but the simple logarithmic 
behaviour Coo ~ ln(l/i) is substituted by the double 
logarithmic singularity 



Coo ~ln l-Fflln 



(11) 



l + ^ln( i 



(9) 



The behaviour of the critical spin-spin correlation func- 
tion is instead predicted not to change in presence of im- 
purities, thus its anomalous dimension rj retains its pure 
Ising value rj — I /4. 

Recent numerical works^^^^^ have provided evidence in 
contrast to what has just been presented above. These 
findings show quantities such as the susceptibility and 
the correlation length to display simple power-law sin- 
gularities, at criticality, with exponents 7 and v varying 
continuously with disorder in a way, however, that their 
ratio "//h' is kept constant at the pure system value. Ac- 
cording to this weak universality^^ scenario, the ratio of 
exponents should not depend on disorder, and the specific 
heat was observed to saturate at a nondivergent value as 
t^O. 

In a finite-size numerical investigation like ours it is 
also essential to translate the predictions of these scenar- 
ios into their finite-size counterparts. This is straight- 
forward and well known in the case of power law singu- 
larities (and power law corrections to scaling), but less 
trivial within the logarithmic corrections scenario. These 
questions are dealt with in detail in the Appendix and, 
where they arise, in the following Section. 

As already stressed in the Introduction, we have tried 
to look at our results with open mindedness: none of the 
two scenarios is chosen as a reference, we rather attempt 
to let our analysis to decide for those of the two which 
consistently fits the whole series of our data. 



IV. MAIN RESULTS 
A. Determination of the critical temperature 

The 2d spin diluted Ising model exhibits a phase 

transition from a paramagnetic to a ferromagnetic low- 
temperature phase, provided the spin density p is larger 
than the critical site percolation density pc — 0.592745.^^ 
The model unfortunately lacks the possibility, known for 
its random-bond counterpart, to exactly determine the 
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critical temperature as a function of the degree of dis- 
order from duality. This is one of the reasons which 
render the bond disordered version more attractive for 
a numerical study, i.e. the possibility of studying its crit- 
ical behaviour by sitting exactly at T^- A very precise 
numerical determination of Tc is thus called for in our 
case, and it has indeed been reached by a joint analysis 
of the finite-size behaviour of the correlation length (as 
defined in Eq. (8)) and of the dom,ain wall free energy 
(or interfacial tension), the latter obtained by compar- 
ing free energies of systems with periodic and antiperi- 
odic boundary conditions in the direction perpendicular 
to the strip. Antiperiodic boundary conditions in a fer- 
romagnetic system force an interface along the cylinder's 
length. In practice, the interface is created by introduc- 
ing an antiferromagnetic seam along the cylinder, i.e. by 
reversing in each row the bond J ^ —J between two 
fixed spins (e.g. spins L and 1). The interface (domain 
wall) free energy per unit length cr^ is given by the dif- 
ference in free energy between the system with periodic 
and antiperiodic boundary conditions, and reads 



0.29 



paL = - In 



71 



(12) 



with ji^abc the largest eigenvalue of the TM with an- 
tiperiodic boundary conditions (compare Eq. (12) to Eq. 
(8))- 

Phenomenological renormalization^^'^'^ predicts that 
the correlation length scales as L at criticality. The 
correct formula reads 



^l' = L-^{A + B^Ly-^ + ...) 



(13) 



where the corrcction-to-scaling terms (the exponent y^^r 
is related to the presence of irrelevant scaling fields and 
it is known to be yirr = —2 for the Ising model^^) split 
the exact crossing of the curves i/^L versus /3J at the 
critical point for different values of L into a sequence 
of distinct intersections of abscissa {PJ)l (see Fig. 1). 
Extrapolating this sequence to its asymptotic value pro- 
vides an accurate determination of the critical tempera- 
ture. The task is accomplished by a combination of dif- 
ferent extrapolation techniques. They are basically: the 
Bulirsch and Stoer algorithm discussed by Henkel and 
Schiitz;^^ the three-point iterated fit method presented 
by Blote and Nienhuis;^'' a third method, close in spirit 
to the latter, which uses fitting procedures at all stages 
of extrapolation.''^ These three different algorithms give 
accurate and consistent results, since they lead to extrap- 
olated values which do not significantly differ from each 
other. They are also used in the following wherever the 
L oo limit of a sequence needs to be extracted. Their 
comparison sets the error bars of our analysis. 

The same scaling argument presented for can be 
applied also to the wall free energy (3a l ~ L~^, and the 
corresponding sequence {l3J)'^ similarly analyzed. 

The apparent "mirror symmetry" of the two sets of 
curves displayed in Fig. 1 suggests also a new method 
to extract f3cJ: the analysis of the sequence (pjy^* of 
mutual intersection points of the curves L/^l and L(3aL 




0.445 



FIG. 1. Set of curves L/n^L (negative slope) and LPol/t^ 
(positive slope) as a function of /3J for different values of L, 
for the pure Ising model. Their mutual intersections, at the 
abscissa {PJYl* , are indicated by open circles, the crossings 
of two Lj-K^h curves for couples of values L, L-\- \ ((/3J)^) are 
drawn as open diamonds, while black dots denote the points 
used in the temperature scan. Is is evident how the former 
sequence is much more rapidly converging than the latter to 
the exact value /3c J = ln(V2 -h l)/2 = 0.440687... . The 
value of the limit of both sequences of points on the vertical 
axis is the exponent rj = 1/4. 



versus /3J for fixed values of L. This sequence also con- 
verges to i3aJ in the limit of large L, but much faster than 
the previous sequences (/3J)i and {PJ)1, and its extrap- 
olation through the methods described above furnishes a 
much more precise determination of the critical temper- 
ature (upper half of Fig. 2). Already for the rather small 
value L = 6, the pure 2d Ising model value of {(5jy^* 
differs from the exact critical value only after the 11th 
decimal digit. Neither this striking result nor the analysis 
of the (/3J)™* sequence are, to our knowledge, present in 
the literature, and we regard them as remarkable findings 
on their own. 

The whole machinery can be carried over to the di- 
luted case. It is worth remembering that in our in- 
vestigation the spin density p is treated as a parame- 
ter fixed to definite values in the range Pc ^ P < 1 to 
scan the phase diagram. The values of p analyzed are 
p = 0.95, 8/9, 0.8, 0.75, 2/3. 

For p ^1 the "mirror symmetry" of Fig. 1 is lost when 
the value of p substantially departs from the pure system 
value. Logarithmic terms would also appear in Eqs. (13) 
in case of the logarithmic corrections scenario^^ (see the 
Appendix and the discussion in Sec. IV B below). In spite 
of that, the sequences (/3J)^'* still converge quite fast to 
their asymptotic value, again allowing extremely accu- 
rate determination of {j3cJ){p) (see lower half of Fig. 2) 
also for the lower values of p, where closeness to the per- 
colation threshold induces stronger finite-size deviations 
from the asymptotic behaviour of the quantities under 
investigation. 

Fig. 3 presents the phase diagram of the 2d site diluted 
Ising model, showing the values of Tc{p)/J obtained in 
this way for strip-widths up to L = 13. 

The critical slope of the curve Tc{p) versus p at p = 1 is 
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FIG. 2. The sequences of finite-size approximations to the 
critical inverse temperature (/3J)^ (circles), (/3J)2 (squares) 
and (/3J)if* (diamonds) as a function of strip- width L for 
p = 1 (the pure Ising case) and for spin density p = 0.75. It is 
apparent that in both cases the last sequence converges faster 
than the other two to the critical inverse temperature of the 
infinite system represented by the dashed line. Here, as well 
as in the following pictures, no error bars are shown if smaller 
than the symbol's size. 



exactly known^s to be Sc = T-^dT^/ dp\p=i = 2/[ln(l + 
+ ^/2/'K)] = 1.564785 . . . ; from an additional cal- 
culation limited to size L = 11 at a value of p very close 
to 1 and from a numerical evaluation of the derivative, 
we obtain an estimate which differs from the exact result 
by approximately 0.01%. 

Let us finally point out that the method just described 
provides a very precise determination of the critical tem- 
peratures for each of our approxirnaMng systems. We 
make no claim, however, that these values coincide with 
those of the fully quenched system: for fixed p they in 
fact turn out to be slightly different from system (a) to 
(dp^ and need not be the same as Monte Carlo data.^** 
The reason is that nonuniversal quantities, such as indeed 
the critical temperature or the value of the percolation 
threshold pc, do depend on the approximating system 
chosen, while universal quantities, such as critical expo- 
nents, central charge, etc. . . , do not. Our confidence in 
this statement, and consequently in the validity of our 
method, will be addressed in Sec. V. 



B. The exponent -q and the ratio j/v 

Having determined with great precision the critical 
temperature of our model, we now turn to the evalua- 
tion of critical exponents. 

The theory of conformal invariance"**^ relates the ampli- 
tude A of the leading term in the FSS behaviour of ^ (Eq. 
(13)) to the anomalous dimension of the spin-spin corre- 
lation function via A — nrj, providing a simple method 
to extract rj just by extrapolating the sequence L/^l to 
L ^ oo. The same procedure can be applied as well to 
L(3aL since the amplitudes for the correlation length and 
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FIG. 3. The phase diagram of the 2d site diluted Ising 
model, with its ferro- and paramagnetic phases. The dia- 
mond and the circle denote respectively the exact percolation 
threshold and the transition temperature of the pure Ising 
model, the squares correspond to the values of spin density p 
analyzed in this investigation. The line is nothing but a guide 
to the eye. 

for the wall free energy are equal. 

In the pure Ising case, 77 is known to be 1/4. Our 
extrapolations unambiguously give the same value 



V = 



(14) 



for all the values of p considered, with a maximum es- 
timated error of 0.0004 (the error arising from the un- 
certainty in the location of the critical temperature be- 
ing much smaller than the difference between final values 
from different extrapolation procedures). From this we 
can conclude that the exponent r] does not depend on dis- 
order, a result which has been known already^'^'^"^'^'*"'^'"^"'^ 
but without the precision shown above. 

Note that the correlation length ^ extracted from the 
ratio of TM eigenvalues (Eq. (8)) within the EE A refers 
to the decay of the average spin-spin correlation function 
in the disordered system. This is to be contrasted with 
the correlation length that describes the typical decay 
of the correlation function in a given disorder configu- 
ration, a quantity readily obtainable from a correspond- 
ing ratio of Lyapunov exponents within a random TM 
approach. ^^'^^ Our identification may be confirmed by 
plotting L/£_L versus l/i^ (Fig. 4) and noticing that no 
corrections to a linear dependence smoothly extrapolat- 
ing to 1 /4 appear, as it could instead be expected for the 
behaviour of the typical correlation length (see in partic- 
ular the discussion focusing on Fig. 1 of Ref. 41, and Ref. 
11). 

We would like to point out that the smallness of the 
estimated error of the exponent 77 is of particular rele- 
vance. It can in fact be also read as a strong confirma- 
tion that our approximations are sufficient to place the 
model under study in the right universality class, and 
that universal quantities therefore turn out to be cor- 
rectly reproduced. 

This is confirmed by the behaviour of the magnetic 
susceptibility data obtained via the group theoretical re- 
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FIG. 4. The scaled inverse correlation length L/tt^l plot- 
ted as a function of for the pure Ising case (circles), for 
p = 0.80 (diamonds) and for p = 2/3 (squares). The lines 
just connect the displayed points. Note the vertical scale. 

duction technique. The finite size scaling analysis of the 

susceptibility x predicts a loading L dependence of the 
form XL ~ L'^/" , and the ratio ^/v is extracted by ana- 
lyzing the sequence 

ln(xi,+i/Xi) 



L In (i + 1/L) 



(15) 



Corrections to scaling generally modify the pure power 
law, and an extrapolation to L ^ oo is needed. Cor- 
rections may be again of pure power law form in L (in 
the conventional situation arising from irrelevant scaling 
fields), or indeed by terms containing logarithms of L (in 
the log-corrections scenario). In the latter case, if the 
largest available system size is too small, the data might 
still reflect a preasymptotic regime, and the very idea of 
extrapolating a sequence like (15) may lose its validity in 
that case. We refer the reader to the Appendix for a dis- 
cussion of this point. We are however not able to detect 
the non-universal corrections predicted to be present in 
the preasymptotic regime of the log-corrections scenario. 
They are either absent, or too small to be detected (Eq. 
A16d). The same is true for 77. 

Extrapolating the sequence (15) provides the result 



7 

V 



(16) 



for all the points studied in the present paper. The es- 
timated error is around 0.001, bigger than for i) because 
a shorter sequence (Lmax = 9) has been utilized. These 
data show that the ratio is independent of the de- 
gree of dilution, a finding in complete accordance with 
the constancy of rj via the Fisher's relation ^/v = 2 — r], 
which is thereby demonstrated to be satisfied with a pre- 
cision of the order of a tenth of a per cent. 

The observed constancy of r; and 7/z/ is unfortunately 
not suflicient to discriminate between the two admissible 
scenarios, since both predict that the ratio retains the 
pure Ising value 7/4, irrespective of disorder. It only 
shows that for these two quantities extrapolations are 
not noticeably misguided by logarithmic corrections, in 
case they are present. 



C. The central chcirge 

Conformal invariance^^ provides a link between the 
values of the dimensionless free energy per site at crit- 
icality /l in a finite strip with periodic boundary con- 
ditions, and the central charge c which characterizes the 
universality class of a conformally invariant model. The 
relation reads 



fL 



f - — + 



(17) 



with /(X) the nonuniversal bulk (infinite system) free en- 
ergy per site. Eq. (17) is supposed to retain its validity 
also in a random system, provided c is substituted by an 
effective central charge c' (Rcf. 33). The ncxt-to- leading 
term, omitted in (17), differs depending on whether log- 
arithmic corrections are present or not. In the latter case 
it is believed to depend on L as B f L~'^'^y^''^ = BfL~^ if 
we adopt the value Uirr = —2 of Ref. 37 or, cquivalcntly, 
that born out from the analysis of Sec. IV B above. In the 
former case, instead, Ludwig and Cardy'^'^ derived the re- 
sults presented in Eqs. (A18) of the Appendix. To avoid 
dealing with the /oo term, we prefer to extract finite-size 
estimates of c' through the formula 



c'l = - {Il+i 



2L + 1 



(18) 



and analyze them both by extrapolating the above se- 
quence to c', and possibly looking for finite-size correc- 
tions (next-to- leading terms). 

The outcome of our analysis does not sensibly differ 
from what we obtained in the study of the exponent 77. 
That is, while c' appears to be clearly given by 



1 



c = 



(19) 



equal to the pure Ising value, for all the points studied 
(with a maximum error of 0.0003), we are not able to in- 
dividuate any correction of the kind of those expressed by 
Eq. (A18). The plots of versus are lines showing 
no appreciable bending. 

The estimates c'j^ converge to their asymptotic value 
always from above, both in the pure as well as in the 
diluted case. This is at least a necessary condition for 
the reflection positivity property, which holds for unitary 
models but it might not hold for random models. The 
discussion on this point in Ref. 43 is indeed based on Eqs. 
(A18) which would predict a convergence from below, 
the flrst correction terms to the value c' = 1/2 being 
negative. But again, they are either absent, or too small 
to be discernible within our errors, or convergence from 
below sets in for system sizes beyond those accessible by 
our TM calculations. 



D. The analysis of the correlation length 

The correlation length provides us with another tool to 
push our investigation further: the study of the exponent 
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TABLE I. Various exponents for different values of spin 
density p. In the second column the exponent v as extrap- 
olated from Eq. (22). The third column shows the value 
of the ratio a/v extracted from Eq. (25). The last column 
provides a test for the validity of the hyperscaling relation 
2/v - ajv = d = 2. 



p 


V 


a 
1/ 


2 a 

V V 


0.95 


1.054 ± 0.002 


-0.102 ± 0.008 


1.999 ± 0.012 


8/9 


1.113 ± 0.001 


-0.214 ± 0.007 


2.011 ± 0.009 


0.80 


1.15 ± 0.03 


-0.30 ± 0.02 


2.04 ± 0.06 


0.75 


1.18 ± 0.02 


-0.35 ± 0.02 


2.05 ± 0.05 


2/3 


1.23 ± 0.05 


-0.445 ± 0.005 


2.07 ± 0.07 



I/, extracted, in addition to 77, from a FSS analysis of ^, 
can be of particular relevance since it is a single isolated 
exponent, therefore it is predicted to show a different 
behaviour in the two scenarios. 

Our analysis follows general reasoning. A recent ap- 
plication to disordered Ising models has been provided 
by Aarao Reis eA, al.}"^ Two quantities at our disposal 
can be subjected to the same investigation, the correla- 
tion length ^ as well as the domain free energy /3cr: both 
show the same asymptotic behaviour, though with differ- 
ent corrections to scaling. 

First we define the derivatives of these quantities with 
respect to temperature: 



m(*) = 



dm 

dt ' 
dt 



(20) 
(21) 



A pure power-law divergence of ^{t) at criticality like 
t~'^ implies ^(t) ^ t""^^, which in turn translates into 
the finite-size prediction ^ L^^^^^ . The same rela- 
tions hold respectively for (3a{t), ii'{t) and /x^. From the 
liL sequence (and similarly form the /U^ sequence), one 
obtains finite-size approximations ul of the correlation 
length exponent via 



In (fiL-L/fJ-L) 

ln{L + l/L) 



(22) 



Extrapolated values clearly appear to vary with p, and 
to change continuously from the Ising value p = 1 to 
the percolation value v = 4/3,^^ in a way shown in Ta- 
ble I. The large error bars on the values correspond- 
ing to stronger dilution are due to degradation in the 
quality of the fits, and even to the appearance of non- 
monotonicities in the sequence i^l versus L extracted 
from the domain wall free energy data, which prevent 
an extrapolation as precise as those for larger p values. 
The cause is probably the closeness of the percolation 
threshold, which requires a corresponding enlargement in 
system size in order to disentangle the system's correct 
critical behaviour from crossover effects. 

A similar data trend was extracted also from Monte 
Carlo simulations,^" and analogous variation with the 
strength of bond disorder was observed in Ref. 12. The 
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FIG. 5. The quantities {pl/L^Y and (^'^/L^)^ at critical- 
ity plotted as a function of In L, for various values of p. 



variation of as a function of p at constant 77 and con- 
stant 7/z/ would qualify the observed non- universality — 
if it must be upheld — to be of the weak form. 

On the other hand, one should be wary of the fact that 
preasymptotic logarithmic corrections might invalidate 
the very idea of extrapolating a vl sequence of limited 
length; sec (Eq. Al6c). This was the attitude advocated 
by Aarao Rcis et al., who did not give much weight to 
an analogous variation of their extrapolated exponent v 
as a function of disorder. Rather they proposed to fit 
the L dependent p^ data to a form expected to hold in 
an intermediate (preasymptotic) range of L values within 
the logarithmic corrections scenario, viz. 



(23) 



which predicts (/xl/L^)^ to be linear in InL. Note that 
this expression used by Aarao Rcis et al. was derived 
from heuristic considerations. It formally agrees with the 
true expression (A16b) only to first order in InL, but it 
predicts the wrong sign for the coefficient of the leading 
In L behaviour of the effective size-dependent correlation 
length exponent ul (Eq. A16c). 

In Fig. 5 we plot our results in the same form as Aarao 
Rcis et al. do, both for p^ and for p'j^. The data show ex- 
actly the same trend as observed in Ref. 12, both with re- 
spect to system size L and with respect to the strength of 
the disorder. From this observation we derive additional 
strong confidence in the validity of our method. The 
pure system behaviour is soon substituted by an appar- 
ently linear term in InL for (/xl/L^)^ in the increasingly 
disordered system, as would be expected from (A16b). 
Notice, once more, the same asymptotic behaviour oi pL 
and p']^, with different corrections to scaling. 

The data can, however, hQ consistently interpreted 
within a scenario of nominivcrsal critical exponents as 
well. Indc;(;d, for a power law divergence of the correlation 
length with v > I, FSS would predict {ph/L'^f ~ L''^ 
with u! = 2 — Ijv at sufficiently large L for the data 
plotted in Fig. 5. In the disorder range considered, lo 
is a rather small quantity as the v extrapolated from 
Eq. (22) only slightly exceed 1. The data might there- 
fore easily look like exhibiting a (1 — A InL) behaviour 
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for the accessible range of system sizes, which was the 
form assumed in Ref. 12. This might however just be the 
result of an expansion to first order in the small quan- 
tity LvlnL {L~'^ ~ 1 — w InL). Upon closer inspection 
a curvature compatible with an L~'^ behaviour is indeed 
discernible in the figure for larger disorder. Using a fit 
against L"" we have another way of analyzing our data 
to determine i/, and the values obtained this way are com- 
patible with those obtained from extrapolation of the 
data. Nonetheless, a curvature of the (/Ul/L^)^ data may 
occur also in the log-corrections scenario (depending on 
the value of Di in Eq. (A13)). 

We can however affirm with certainty that neither 
Aarao Reis et al. nor we are seeing preasymptotic ef- 
fects of logarithmic corrections in the data for vl, be- 
cause they are increasing with L (in our case at least 
for p > 0.8), whereas according to FSS they should de- 
crease; see Eq. (A16c). Thus, if logarithmic corrections 
are present, they are for the available system sizes still 
masked by other corrections to scaling. Again, it seems 
that further investigation is required in order to be able 
to take sides on the critical behaviour of this system. 

E. The specific heat 

The logarithmic correction scenario predicts for the 
specific heat the double logarithmic divergence of Eq. 
(11). This translates, at criticality, into the FSS 
prediction^ 

~C^ + C; In (1 + InL) , (24) 

where the pure system critical behaviour is recovered by 
the vanishing of C2 (but with C[C2 = const.) for p = 1. 
On the other hand, a scenario of nonuniversal critical 
exponents would lead to 

Cl^Coo+ C^i"/" . (25) 

Though a/u is a ratio of exponents, it is expected to vary 
with p, being directly related to v via the hyperscaling 
relation 2 — a = du, with d the system's dimensionality. 
Moreover, this relation imposes a severe constraint on a: 
since v is found to be greater than 1, a should be smaller 
than 0, in other words the specific heat should then turn 
out to be nondivergent, with its finite-size estimates sat- 
urating to a value Coo at 

We plot our data for the specific heat in the same 
form as in Ref. 13, i.e. against both InL and InlnL 
(Fig. 6). The pure system divergence is well reproduced 
by the straight line the data show when plotted against 
In L, equivalently by the upward curvature when against 
InlnL. Exactly as it happens in Ref. 13, as disorder is 
switched on and increased, it is apparent that both curves 
tend to bend downwards, the former markedly deviating 
from, the latter instead approaching, a straight line (from 
top left to bottom right of Fig. 6). The law describing 
these finite size data has therefore to be at least less di- 
vergent with L than InL. However, the same kind of 
behaviour must at moderate system sizes be expected in 
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FIG. 6. The finite-size approximants of the specific heat at 
criticality, Cl , versus In In L (circles) and In L (squares) . 

both scenarios (as we shall argue in greater detail in the 
concluding section). It is therefore difficult to conclude 

definitely in favour of Eq. (24) : here more than elsewhere 
the need for significantly larger system sizes is particu- 
larly felt. 

If we fit the data to a nondiverging size dependence, 
in accordance with Eq. (25), there is a consistency check 
to carry out, i.e. one may ask oneself whether the hy- 
perscaling relation is verified or not, given the f values 
from the correlation length. The fitting of the specific 
heat to the form (25) is performed first by taking all the 
points into account, then by reducing the data set by 
successively increasing the initial value of L considered. 
In this way one can both check for the stability of such 
a fit and get a sequence of values of a/v which may be 
subjected to subsequent extrapolation procedures, when 
needed. The rcsiilts of this analysis arc shown in Table 
I. Its main outcome is the fact that, within a weak uni- 
versality interpretation of the data, a small, seemingly 
systematic offset from hyperscaling, though still within 
estimated error bars, appears. 

We have tried to discriminate between Eqs. (24) and 
(25) on the basis of a test. Although such a test ap- 
pears to favour the log-corrections scenario by roughly 
an order of magnitude in x^, the values for both al- 
ternatives are so small (typically 0(10^^) and 0(10"'^) 
respectively for a choice of the data set) as to make it 
doubtful whether a meaningful model selection should be 
based on them. In particular, our data are not normally 
distributed random data (indeed they are deterministic), 
as instead strictly required by a reliable analysis. 

V. DISCUSSION AND CONCLUSIONS 

In this paper we studied the critical behaviour of the 
2d Ising model with quenched random site dilution via 
the EEA. Our purpose was twofold: first we wanted to 
address again the question about the reliability of the 
method; second, we extended the study started in Ref. 21 
to systems of larger sizes, also in view of the new results 
which have appeared in the literature since then,^^~^^ in 
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order to allow perhaps a clearer discrimination between 
the two contradicting pictures of the system's critical be- 
haviour that survived so far, the logarithmic corrections 
and the weak universality scenario. 



A. Reliability of the Method 

As the description of quenched disorder within our 
EE A is only approximate, we have to worry about how 
good it actually is. Before gathering the different pieces 
of evidence accumulated during our numerical study con- 
firming the validity of our approach, we would like to 
address the question from a more general viewpoint. In 
Ref. 21 the correlation length ^j, describing the asymp- 
totic decay of the correlation function {rij ) of the dis- 
order degrees of freedom was introduced and studied. Its 
finite-size estimates can be computed from the ratio of 
two TM eigenvalues in complete analogy to (8) by replac- 
ing 72 with 72 , where 72 is the second eigenvalue of the 
symmetric block of the transfer matrix F. 

As regards the analysis of this new correlation length, 
we have both bad and good news. The bad news is that 
^fc actually diverges at criticality in our different approx- 
imating systems. This escaped our attention in Ref. 21: 
while was observed never to exceed a few lattice spac- 
ings for the accessible system size, its size dependence 
at Tc had not been monitored. However, before con- 
cluding from this that our systems provide only a rather 
poor description of quenched disorder, we have also deter- 
mined the behaviour of the ratio Rl = £,k,L{Tc)/^L{Tc) 
of these two correlation lengths at criticality, and the 
results of this study may be taken for the good news. 
We find that Rl 1/8 at large L, independently of p. 
Since this limit represents nothing but the ratio rj/rik of 
the anomalous dimensions of the spin and occupation- 
variable correlation functions at Tc, this implies that 
rjk = 2, hence unusually large. The fc^-correlator is thus 
'almost summable' at criticality (its sum has a logarith- 
mic infrared divergence). Note that the ratio Gk{r)/G{r) 
behaves as 

Gk{r) J_ 



G(r) 



(26) 



and thus decays to zero, stating that ki correlations are 
negligibly small if compared to crj correlations at large 
distances. This is why the systc^m may be regarded as 
effectively quenched despite remaining correlations be- 
tween the ki, and it may be regarded as the reason why 
our results compare so favourably with those obtained 
via more conventional approaches. 

Indeed, both checks against exact results, wherever ob- 
tainable, and direct comparison of our data with those 
obtained by the random TM approach (which allows for 
an in principle exact treatment of the disorder) increase 
our confidence in the validity of our method beyond any 
reasonable doubt. Among the exact results correctly re- 
produced by our method we mention: i) data obtained 
for a one-dimensional system;**^ ii) the value of the ini- 
tial slope of the critical line in the phase diagram, re- 
produced to within 0.01%; in) the correct value of the 



connectivity length exponent Vp and of the crossover ex- 
ponent 4> = ^p/i^ = 1 at percolation;^^ v) the precision 
with which the values of the exponent 77, of the ratio 
7/1/, and of the central charge are determined for all 
the points investigated in the phase diagram. In ad- 
dition, as mentioned above, the random TM data (ob- 
tained in the case of bond disordered systems)^^'^^ ex- 
hibit the same finite-size signature, and qualitatively the 
same behaviour as regards their dependence on the dis- 
order strength as those obtained in the present paper. 

A definite advantage of our method over other numer- 
ical methods dealing with disordered systems, is that it 
does not suffer from non-sclf-averaging difficulties, and 
provides directly through the ratio of TM eigenvalues the 
average correlation length, not the typical one, such that 
the successive analysis turns out to be more straightfor- 
ward than, e.g., in Ref. 12. 

The results discussed so far were all obtained within 
approximating system (d). Comparative studies includ- 
ing results also from the other systems (a) -(c) were pre- 
sented in Ref. 21. They showed that systems (a)-(d) 
appear to be in the same universality class as regards 
their critical behaviour. Non-universal quantities, how- 
ever, such as the critical temperature itself, or the value 
of the percolation threshold Pc, do depend on the approx- 
imating system, differing between those with and without 
the plaquette constraint and those estimated by other nu- 
merical works. ^^'"^^ The conclusion is that the description 
of quenched disorder, though only approximate, appears 
to be precise enough to put all our approximating systems 
into the same universality class as the the fully quenched 
system, and this is all that is required for the purpose of 
the present study. 



B. Discriminating Between the Scencirios 

We now come to our second point, concerning the dis- 
crimination between the two scenarios proposed for the 
critical behaviour of the model considered in this paper. 
Our view is that some of the FSS investigations which 
have appeared in the literature - including the recent 
ones^'^°'^^~^*'^^'^°'^^ - may perhaps not be as decisive on 
this matter as their authors have tended to believe. We 
shall now discuss some issues on which, wc believe, clar- 
ification or amendments are needed, presenting a critical 
review of available material in the light of our own data 
analysis. 

The results of Refs. 20 and 21 have been interpreted as 
providing evidence in favour of a weak universality sce- 
nario. The former is a Monte Carlo investigation both at, 
and off criticality of the site diluted system. On the basis 
of data analyses, Kim and Patrascioiu concluded that 
their off-critical simulations of susceptibility and correla- 
tion length were better described by modified power laws 
than in terms of the log-corrections scenario. No analo- 
gous discrimination was attempted for the specific heat, 
and hence no check of hyperscaling was performed. The 
main weakness of their off-critical simulations is that re- 
duced temperatures are still sizeable, and the constancy 
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of the ratio 7/1^ appears slightly less well satisfied than 
the individual errors on the values of 7 and v would 
allow. Simulations at criticality yielded a specific heat 
very slowly increasing with system size, and arguments 
in favour of a saturation were advanced. The constancy 
of 7/1^ was clearly shown and the validity of the Fisher 
relation confirmed, but - as noted above - this cannot be 
taken to support either scenario over the other. 

Our own earlier TM strip-scaling results, ^""^ too, were 
interpreted in terms of weak universality. Possible effects 
of logarithmic corrections in the effective size-dependent 
exponents vj^ were looked for, but were either absent or 
simply not discernible due to the rather moderate strip 
widths available in that study. Indeed, the rather limited 
system sizes may be taken as one of the weakest points 
of that investigation. Moreover, specific heats had not 
been computed and so hyperscaling within a varying ex- 
ponents picture was not checked. Again, the observed 
constancy oi^jv and may bo regarded as a piece of ev- 
idence for the reliability of the method, but not in favour 
of either scenario, except in so far as implying, that z/ 
critical behaviour in the model were non-universal, the 
observed non-universality would have to be of the weak 
form. 

Let us now turn to the recent investigations which have 
been taken to support the log corrections scenario. They 
are either of the TM strip-scaling^°'^^'^^ or of the Monte 
Carlo«'i4.i6 type. 

The strip scaling data of Aarao Reis et aL^°'^^ for 
the correlation length in the bond disordered system ap- 
peared to give values of v slightly greater than that of the 
pure system, which the authors, however, discarded, at- 
tributing them to preasymptotic effects originating from 
logarithmic corrections. Their interpretation of the data 
in these terms does, indeed, describe the data rather well 
on the level of the /Ul, but not on the level of the size 
dependent exponents vj^. Moreover, we recall our discus- 
sion of this point in Sec. IV D, according to which the 
power law picture provides a consistent interpretation of 
the data as well. 

Turning to the critical specific heat data presented in 
Refs. 12 and 13, they were interpreted as "clearly sug- 
gesting a divergence in the thermodynamic limit" . This 
point was claimed to be strengthened by pushing sys- 
tem size up to L = 18 (and in less precise simulations 
even to impressive L = 23), and by noting that spe- 
cific heat data appeared to be halfway between single- 
(InL) and double- (InlnL) logarithmic dependence on 
strip width, from which a divergence in the thermody- 
namic limit was inferred. One should, however, note that 
at moderate L this kind of behaviour is to be expected in 
hoth scenarios if the correlation length exponent is close 
to 1 and, hence, a only slightly negative. The difference 
is that Cl will continue to increase at least as fast as 
Inlni for all L, in the log-corrections scenario, whereas 
there will be a crossover to a growth which is slower than 
InlnL at L, ~ exp{— i^/a} in case of a modified power 
law (25) with a < 0. Tentatively accepting the value 
V ~ 1.083 determined by Aarao Reis et al}° for the 
r = J1/J2 = 0.25 case investigated in Ref. 13, one would 



have to locate the crossover length at ~ 680. Thus 
Lmax = 23 is still much too small to allow concluding in 
favour of either scenario. 

We have carried out data analyses on the raw data 
of Ref. 13, comparing the two scenarios, as we did on 
our own data. Power law fits and fits according to (24) 
show no significantly different quality. In particular, it 
seems that a double logarithmic law provides a better 
fit to the data selected from a window [Lmin,23] with 
-^min up to 8, while power laws give a smaller for 
the larger imin up to 13; for still larger L the result of 
the fits become questionable for both hypotheses. Other 
kinds of fits have been tried, e.g. by selecting a movable 
window [Lmin, ^min + 7] and letting Z/min run over the 
data, but without any significant improvement. Only by 
discarding the larger L values, which are still much too 
noisy (as it is obvious by a quick look at the derivative 
information plotted in Fig. 3 of Ref. 13), one gets values 
of slightly in favour of the double logarithmic form 
(24). Incidentally, however, the power- law fits lead to an 
estimate for a/v compatible via hyperscaling with the v 
value reported in Ref. 10. This analysis cannot thus be 
seen as conclusive. 

The Monte Carlo study of site diluted system by 
Ballestcros et al}^ shows critical specific heat data rea- 
sonably well fitted (in terms of y^) by a double logarith- 
mic form, at least for intermediate dilution. The same 
is true for a recent study of Selke et al}^ However, no 
power-law fitting is attempted for comparison. In earlier 
work on bond disordered systems,^ such power-law fits 
had been attempted, and were regarded inferior to dou- 
ble logarithmic ones. Note, however, that the maximum 
system sizes studied in Ref. 7 are for r = 0.25 still below, 
and for r = 0.1 above but still very close on a double log- 
arithmic scale to the crossover lengths expected from 
the V values reported in Ref. 10. 

The vl data presented in Ref. 14 have, in our view, 
somewhat too large error bars to allow concluding with 
confidence that vl ^ ^ for large L for all the values of p 
also because Ballesteros et al. force the intercept through 
1 rather than fitting it. The reasonably largo values of 
L reached in this investigation would probably set the 
system in the asymptotic regime (at least for the stronger 
disorder), and the FSS expression used by Ballesteros et 
al. is indeed = 1 + ^'/lnL. Nonetheless the constant 
A' in this expression should, strictly speaking, come out 
to be 1/2, independently of the degree of disorder (see 
Eq. A17c), a property which their data and their fits do 
not respect. 

In the off-critical Monte Carlo study of Talapov and 
Shchur^, bond disorder was observed to lead to increased 
values for the magnetization and susceptibility exponents 
f3 and 7, while the behaviour of the specific heat showed 
good agreement with the double logarithmic form (11). 
Together with the Rushbrooke equality a + 20 + 'y = 2, 
these findings embody an inconsistency from which the 
authors concluded that their increased values for /3 and 
7 cannot be asymptotic. On the other hand, it turns out 
that the specific heat data of Ref. 8 can be fitted equally 
well to a power law form with a value of a compatible 
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via the Rushbrookc relation with the increased values of 
P and 7.''® Note also that the ratio P/j in the disordered 
system is the same (to within a fraction of a percent) as 
in the pure system, as would be expected in a weak uni- 
versality scenario, whereas individually (3 and 7 change 
in the 6% range. 

The difficulties of FSS are avoided in the series ex- 
pansion study of Roder et al.;^^ their determination of 
the exponent 7 in Eq. (10), which appears to saturate at 
7/8 for sufficiently strong disorder, may perhaps be taken 
as the strongest piece of evidence currently available in 
favour of the validity of the logarithmic corrections sce- 
nario. Still, the analysis suffers from the relatively small 
length of the series, at least in some regions of the phase 
diagram, and one would wish this to be made more con- 
clusive by including higher order terms, so as to reduce 
error bars. It is also worth noting at this point that their 
claim that Xz,(T'c) is unaffected by logarithmic corrections 
cannot be upheld: see Eq. A15. 

We now turn to the present investigation. First, by 
combining finite size signatures of correlation length and 
domain wall free energy, we have in particular been able 
to locate critical temperatures with extreme precision. 
Second, we have significantly enlarged our system sizes. 
In interpreting our data, maximum care was constantly 
taken to be open in both directions. Based on results 
of Cardy and Ludwig, a more systematic FSS analysis 
than in previous numerical studies has been performed. 
As we will presently show, neither way of looking at the 
available finite size data is completely satisfactory. 

The value of the effective central charge is foimd to be 
c' = 5 with extremely high precision. If conformal invari- 
ance and reflection positivity would also hold for disor- 
dered systems, this finding would put the model into the 
Ising universality class and exclude continuously varying 
exponents. However, the results of Ludwig and Cardy^^ 
imply that the latter condition does not hold, at least for 
the weakly disordered ferromagnetic Ising model. 

Considering the correlation length exponent ly, we ob- 
serve that our FSS estimates converge to values contin- 
uously varying with p. One might suspect, of course, 
that our extrapolations are misled because the algorithm 
would not pick up slowly varying logarithmic corrections. 
However, if we adopt this hypothesis, we find it not easy 
to reconcile with the fact that not even the slightest such 
effect is detectable in our extrapolations of r/, j/iy and 
the central charge c, even though similar, albeit smaller 
offsets must be expected to occur in these quantities as 
well, as is borne out by our analysis of FSS in the Ap- 
pendix. 

The fact that the critical specific heat data may be 
well fitted both to a double logarithmic behaviour (24) 
and to a cusp singularity (25) shows that it is difficult to 
discriminate with confidence between the two laws. How- 
ever, we recall that a fit according to a modified power 
law (25) entails a weak violation of hyperscaling, though 
still within error bars. 

In conclusion, as regards the two issues raised in this 
paper, we believe to have provided, if not a proof, sat- 
isfying arguments in favour of the trustworthiness and 



accuracy of the EEA to disordered systems. As to the 
second issue, we feel that further investigations would 
still be needed to provide a clearer and definite evidence 
in favour of either picture. The problems seen in the FSS 
analyses of the system are not exclusive to our results, 
but are - as we believe to have demonstrated - common 
to virtually all previous FSS data obtained so far. 

Let us finally emphasize that we have no reasons to a 
priori distrust the correctness of the theoretical picture 
that began to emerge through the work of Dotsenko and 
Dotsenko^ and the improved and corrected versions of 
Shalaev^, Shankar^ and Ludwig^. Still, in an ideal world, 
one would like to see this picture supported by numerical 
evidence much better than currently available. 
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APPENDIX: FSS WITH LOGARITHMIC 
CORRECTIONS 

In the present Appendix we collect the main results 
of the renormalization group (RG) analysis for the case 
where small amounts of disorder constitute a marginally 
irrelevant perturbation at the pure system's fixed point, 
thus giving rise to logarithmic corrections, as in the field 
theoretical approach to the weakly disordered 2d Ising 
model. We do this both for the sake of completeness, 
and to obtain sound understanding of the effects such 
corrections would have in a FSS analysis. This will turn 
out to be relevant not only in the asymptotic regime L 3> 
1 but especially in the preasymptotic regime as specified 
below. 

We feel it particularly important to collect these re- 
sults here, because the literature in the field abounds 
in heuristic derivations which have sometimes produced 
erroneous results and still more frequently even miscon- 
ceptions as to what the effect of logarithmic corrections 
might be in the FSS signature of various thermodynamic 
functions and critical exponents. 

The following is basically exploring results Ludwig and 
Cardy^"^ obtained in a replica approach to the problem, 
which are in turn based on earlier results on logarith- 
mic corrections obtained by Cardy in a more general 
context. 

If g denotes the marginally irrelevant coupling arising 
from disorder averaging in the (replicated) Hamiltonian, 
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and by the set of other scahng fields at the pure 

systems critical point, the RG flow-equations read 



dg 
dl 

dl 



: -nbg^ + 0{g^) (Al) 
: UnUn - 2nbngun + 0{g'^Un) , (A2) 



where the t/„ are the eigenvalues of the linearized RG 
equations and the 6„ arc operator product expansion 
(OPE) coefficients, which can be formulated in terms of 
corresponding three-point functions. 

Up to the order shown, these equations integrate to 



with go = g{0) and 

u„(Z) = ^i„(0)e^"'(l + TTbgol)-^''"/'' = w„(0){*„(0 . (A4) 

Among the {un}, the two relevant fields are Ue which 
couples to the energy density and has RG eigenvalue = 
1, and Ua- coupling to the magnetization density with 
Ucr = 15/8. One may thus put Ug{0) ~ t and Ua{Q) — h. 
In a finite system of linear extent L one may formally 
interpret as another relevant scaling field with RG 
eigenvalue 1. 

In the context of the replica formulation, the OPE co- 
efficients and the coupling g are normalized such that 
the zero-replica limit of the product bgo and of the ra- 
tios bn/b exist and arc finite. The results of Ludwig and 
Cardy imply bgo = 8 A and be/b = 1/4 in this limit, 
whereas b„/b = 0. Here A denotes the (bare) disorder 
strength: A (x p{l — p) for a randomly site diluted model 

The RG equations for the singular part of the free en- 
ergy and the correlation length read 

fit, h, go, L-') = e-'^^f{u,il),u,{l), g{l),o'L-^) (A5) 

and 

i{t,h,go,L-')=B'^{uS),u,{l),g{l),e'L-^) . (A6) 

Irrelevant scaling fields providing additional corrections 
to scaling have for simplicity been suppressed in these 
expressions. As usual, critical behaviour characteristic 
of the infinite system is obtained by considering / and 
^ and their derivatives with respect to temperature (and 
field) at a scale I chosen such that Us{l) = ±1. For \t\ <C 1 
this gives 



(A3) 



1 



1 , '^^90 , 1 
, ■ H In T-, 

t\ \ Ve \t\ 



2be/b 



(A7) 



Inserting the values for the OPE coefficients reported in 
Ref. 33, and i' = 1/ye = 1 for the pure 2d Ising model, 
this gives the divergence of the correlation length at /i = 
to leading order 



m 



87r A In ^ 
t 



1/2 



(A8) 



The fact that ba^/b = entails that magnetization and 
susceptibility will to leading order scale as pure powers 
of the correlation length, in agreement with earlier results 
of Shalaev^. 

Phenomenological renormalization is based on analyz- 
ing the finite size signatures of (A6) and of its tempera- 
ture derivative 

fJ-Lit) = ^ £,{t,h,go,L-'^) 

= e' we(0 Uue{l),u,{l),g{l),e^L-^) , (A9) 

where designates the partial derivative of the 
correlation-length scaling function (A6) with respect to 
Us{l). Finite size signatures at criticality {t = 0, h = 0) 
are obtained by analyzing these quantities at a scale 
e' = Z/, giving 



^l^=L-^ r'(0,0,5(lni),l) 
= i-i$(5(lni)) 



(AlO) 



and 



ML = i^(l + 87rAlnL)-V2^,(0,0,p(lnL),l) 

= i^(l-F87rAlnL)-^/2$e(5(lni)) (All) 

respectively. Eqs. (AlO) and (All) define the univer- 
sal scaling functions <E>(.t) and $e(a;) of the argument x 
which, in the zero-replica limit, is a; = A/(l -|-87rAlnL). 
In the weak disorder limit or at large L this quantity 
is small, and one assumes that an expansion of these 
scaling functions in x exists. For $ one gets ^{x) = 
(jjo -\- (f^ix + (j)2X^ + 0{x^), the coefficients being known to 
be ^0 = '^V from conformal invariance,^° and = on 
account of the vanishing of be,. Thus 

g'°^~'h+^(( i+8ii„J ') 

A similar expansion is expected to hold for $£, i.e., 

^e{x) = (f>eO + <t>elX + 0{x^), SO (with Di = (j)el/(j)eo) 



ML = </'eoi'(l + 87rAlnL)-^/' 



1 + 



LiiA 



+0 



1 + SttA In L 



1 -|-87rAlni 

(A13) 



though the coefficients (j)eO and (f)ei are to the best of our 
knowledge not known. Interest in this quantity derives 
from the fact that, within the phenomenological renor- 
malization group scheme, the correlation length exponent 
y is obtained by computing the sequence of finite size ap- 
proximants vl, defined by 



din 



^ dlnL 



(A14) 



and by extrapolating this sequence to the large system 
limit. Eq. (22) is just the finite diff'erence approximation 
to (A14). 
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Similarly, the FSS analysis for the susceptibility gives 

Xl = L'^''' /<.a(0,0,5(lnL),l) 

= L'^'^^M^L)) (A15) 

with the same coiivcintions for the notation. The ra- 
tio is that of the pmc system, i.e., ^/v = 7/4. 
On account of the vanishing of no additional log- 
arithmic terms above those following from the expan- 
sion of appear in (A15), in contrast to what happens 
for As before, one assumes that an expansion of 
the form = '^xo + 0xi^ + 4>x'2^'^ + 0{x^) exists, 

and defines the effective size dependent ratio (7/2^) l by 
{ll^^L = 91nxi/91nL. It will be useful to introduce 
the abbreviation E-i = (f)^i / (6^0 below. 

For the interpretation of numerical data it is relevant 
to note results for these quantities both in the preasymp- 
totic regime SttA In L -C 1 and in the asymptotic regime 
SttAIuL > 1. 

(i) In the preasymptotic regime, they read 

^-1 =L-i (7r?7-F<^2A2(l- 167rAlnL)) (A16a) 
IJ,L = (peoL'^{l + DiA-{l + 2DiA)4TrAlnL) (A16b) 
UL = 1 + 47rA + (47rA)2 (1 -|- £>i/27r - 2 In i) (A16c) 

O L " i " STrA^iJi + SttA^ (e^ - 2E2 

-|-167r£;ilni) (A16d) 

where we have kept the lowest order in A In X. Additional 
terms down by further factors of A or AlnL are not 

shown. These results exhibit non-universal corrections. 

(ii) In the asymptotic regime, on the other hand, one 
has 



3\ 1 



Ml = ' 



toi^(87rAlnL) 

1 + ^^ + 



InLJ J 



l'L = l + 



21nL 



-I- 



(-) 



7 p 1 
4 '87r(lnL)2 



1 n3 
hTL 



(A17a) 

(A17b) 
(A17c) 
(A17d) 



so the asymptotic corrections to ^j^^, vi, and turn 
out to be universal. 

Ludwig and Cardy have also reported correspond- 
ing FSS expressions for the central charge, in both 
regimes. These are 



c = - - 1287r3 A^ (1 - 247rA In L) (AlSa) 



and 



c=---(lnL)-3 +o(—^ , , 
2 4^^ ISTTlni; ' 



respectively, where we have once more kept the lowest 
order in A In L in the preasymptotic regime. 

Additional terms further down by powers of L^°(l -|- 
STrAlni)"^''"/'', with < —2, will appear due to irrel- 
evant scaling fields. 

Eqs. (A16) - (A18) are of great importance when ana- 
lyzing FSS data, if logarithmic corrections are expected 
to be present. Especially in a strip-scaling approach such 
as ours, the preasymptotic results may turn out of partic- 
ular relevance, since the maximum size available is such 
that one may not reach the asymptotic regime, which is 
likely to be true at least for weak disorder. This has, in- 
deed, been the point of view emphasized by Aarao Reis 
et al.i2 



* E-Mail: mazzeo@pooh.tphys.uni-heidelberg.de 
^ E-Mail: kuehn@tphys.uiii-heidelberg.de 

^ A. B. Harris, J. Phys. C 7, 1671 (1974). 

^ Vik.S. Dotsenko and Vl.S. Dotsenko, Adv. Phys. 32, 129 

(1983). 

^ B. N. Shalaov, Sov. Phys. Solid State 26, 1811 (1984). 

" R. Shankar, Phys. Rev. Lett. 58, 2466 (1987). 

^ A. W. W. Ludwig, Phys. Rev. Lett. 61, 2388 (1988). 

^ B. Derrida, B. W. Southern, and D. Stauffer, J. Phys. 

(Paris) 48, 335 (1987). 
^V. B. Andreichenko, VI. S. Dotsenko, W. Seiko, and J.- 

S. Wang, Nucl. Phys. B 344, 531 (1990); J.-S. Wang, W. 

Selke, and V. B. Andreichenko, Physica A 164, 221 (1990). 

* A. L. Talapov and L. N. Shchur, J. Phys. C 6, 8295 (1994). 

® W. Selke, L. N. Shchur and A. L. Talapov, in Annual Re- 
view of Computational Physics I, edited by D. Stauffer 
(World Scientific, Singapore, 1994), p. 17. 

^° F. D. A. Aarao Reis, S. L. A. de Queiroz, and R. R. dos 

Santos, Phys. Rev. B 54, R9616 (1996). 
" S. L. A. dc Queiroz, J. Phys. A 30, L443 (1997). 

F. D. A. Aarao Reis, S. L. A. de Queiroz, and R. R. dos 

Santos, Phys. Rev. B 56, 6013 (1997). 

D. Stauffer, F. D. A. Aarao Reis, S. L. A. de Queiroz, and 

R. R. dos Santos, Int. J. Mod. Phys. C 8, 1209 (1997). 

H.G. Ballesteros, L.A. Fernandez, V. Martm Mayor, A. 

Mufioz Sudupe, G. Parisi, and J.J. Ruiz-Lorenzo, J. Phys. 

A 30, 8379 (1997). 

A. Roder, J. Adler, and W. Janke, Phys. Rev. Lett. 80, 

4697 (1998); Physica A 265, 28 (1999). 
'-'^ W. Selke, L. N. Shchur and O. A. Vasiliev, Physica A 259, 
388 (1998). 

" V. N. Plechko, Phys. Lett. A 239, 289 (1998). 

R. Kiihn, thesis, Kiel University, 1987. 
" M. Fahnle, T. Holey and J. Eckert, J. Magn. Magn. Mater. 

104-107, 195 (1992). 
^° J.-K. Kim and A. Patrascioiu, Phys. Rev. Lett. 72, 2785 

(1994); Phys. Rev. B 49, 15764 (1994). 

21 R. Kiihn, Phys. Rev. Lett. 73, 2269 (1994). 

22 H. O. Hcucr, Europhys. Lett. 16, 503 (1991). 

A. J. F. de Souza and F. G. Brady Moreira, Europhys. 
Lett. 17, 491 (1992). 
^* K. Ziegler, Nucl. Phys. B 344, 499 (1990); Europhys. Lett. 
14, 415 (1991). 



(AlSb) 



R. B. Stinchcombc, in Phase Transitions and Critical Phe- 
nomena, edited by C. Domb and J. L. Lebowitz (Academic 
Press, London, New York, 1983), Vol. 7, p. 151. 
2^ T. Morita, J. Math. Phys. 5, 1401 (1964). 



CRITICAL BEHAVIOUR OF SPIN DILUTED ISING. 



G. Sobotta and D. Wagner, Z. Phys. B 33, 271 (1979). 

R. Kiihn, Z. Phys. B 100, 231 (1996). 
2^ M. P. Nightingale, Physica A 83, 561 (1976); J. Appl. Phys. 

53 7927, (1982). 
3° M. P. Nightingale, Proc. Kon. Ned. Akad. Wet. B 82, 235 

(1979); 82, 245 (1979); 82, 269 (1979). 

M. Suzuki, Progr. Theor. Phys. 51, 1992 (1974). 

32 J. L. Cardy, J. Phys. A 19, L1093 (1986). 

33 A. W. W. Ludwig and J. L. Cardy, Nucl. Phys. B 285 
[FS19], 687 (1987). 

3^ R. M. Ziff and B. Sapoval, J. Phys. A 19, L1169 (1986). 
3^5 H. W. J. Blote and M. P. M. den Nijs, Phys. Rev. B 37, 
1766 (1988). 

3^ M. Henkel and G. Schiitz, J. Phys. A 21, 2617 (1988). 
3^ H. W. J. Blote and B. Nienhuis, J. Phys. A 22, 1415 (1989). 
3* G. Mazzeo (unpublished). 

3^ M. F. Thorpe and A. R. McGurn, Phys. Rev. B 20, 2142 
(1979). 

^° J. L. Cardy, J. Phys. A 17, L385 (1984). 

*^ S. L. A. de Queiroz, Phys. Rev. E 51, 1030 (1995). 

*2 H. W. J. Blote, J. L. Cardy and M. P. Nightingale, Phys. 

Rev. Lett. 56, 742 (1986); I. Affleck, ihid. 56, 746 (1986). 
*3 J. L. Jacobsen and J. L. Cardy, Nucl. Phys. B 515, 701 

(1998). 

M. P. M. den Nijs, J. Phys. A 12, 1857 (1979). 
*5 R. Kiihn (unpublished). 
"'^ A. Talapov (private communication). 

D. Friedan, S. Qiu and S. Shenker, Phys. Rev. Lett. 52, 



1757 (1984). 



